/** Fichier implémentant un hyperboloïde infini.
 *
 * Ce fichier contient l'implémentation de la classe BaseHyperboloide.
 *
 * @file BaseHyperboloide.cpp
 * @author Johan "Solo" Luisier
 *
 * @date 4 Juillet 2011
 */
#include "BaseHyperboloide.hpp"

#include "Algebre/Solveur.hpp"

#include "Utils/Erreur.hpp"

#include <map>

using namespace std;

/** Constructeur par défaut.
 *
 * Ce constructeur crée un hyperboloïde dont le centre de gravité est à
 * l'origine. Les axes de l'instance sont confondus avec \f$ 0x, 0y, 0z \f$. La
 * section à \f$ z = 0 \f$ est une ellipse.
 */
BaseHyperboloide::BaseHyperboloide()
    : FormeInfinie( "BaseHyperboloïde" ), Sommet( 0., 0., 0. ),
      AxeX( 1., 0., 0. ), AxeY( 0., 1., 0. ), AxeZ( 0., 0., 1. ),
      EchelleX( 1. ), EchelleY( 2. ), EchelleZ( 4. )
{
    ReduitCoordonnes = MatriceGL3( 1. / EchelleX, 1. / EchelleY, 1. / EchelleZ );

    MoteurResolution = new Solveur();
}

/** Constructeur standard.
 *
 * Constructeur permettant d'initialiser tous les champs.
 *
 * @todo Il n'est pas utile de donner l'axe Oz explicitement, il suffirait de
 * donner Ox et Oy.
 *
 * @param[in] sommet centre de l'hyperboloïde.
 * @param[in] axeX axe 0x dans les coordonées locales.
 * @param[in] axeY axe 0y dans les coordonées locales.
 * @param[in] axeZ axe 0z dans les coordonées locales.
 * @param[in] echelleX facteur d'échelle sur l'axe Ox local.
 * @param[in] echelleY facteur d'échelle sur l'axe Oy local.
 * @param[in] echelleZ facteur d'échelle sur l'axe Oz local.
 */
BaseHyperboloide::BaseHyperboloide( const Coord& sommet,
                                    const Vecteur& axeX, const Vecteur& axeY,
                                    const Vecteur& axeZ,
                                    const double& echelleX,
                                    const double& echelleY,
                                    const double& echelleZ )
    : FormeInfinie( "BaseHyperboloïde" ), Sommet( sommet ),
      AxeX( axeX ), AxeY( axeY ), AxeZ( axeZ ),
      EchelleX( echelleX ), EchelleY( echelleY ), EchelleZ( echelleZ )
{
    MatriceSO3 tentative;

    try
    {
        tentative = MatriceSO3( AxeX, AxeY, AxeZ );
    }
    catch( Erreur& e )
    {
        cerr << "BaseHyperboloïde :  Impossible d'initialiser la rotation."
             << endl;
        exit( 1 );
    }

    ReduitCoordonnes = MatriceGL3( 1. / EchelleX, 1. / EchelleY, 1. / EchelleZ );

    VersLocal = Transformation( tentative, -Sommet );
    VersGlobal = VersLocal.inverse();

    MoteurResolution = new Solveur();
}

/** Constructeur de copie, permettant de dupliquer une instance de
 * BaseHyperboloide.
 *
 * @param[in] par instance de Hyperboloide à copier.
 */
BaseHyperboloide::BaseHyperboloide( const BaseHyperboloide& par )
    : FormeInfinie( par ), Sommet( par.sommet() ),
      AxeX( par.axeX() ), AxeY( par.axeY() ), AxeZ( par.axeZ() ),
      EchelleX( par.echelleX() ), EchelleY( par.echelleY() ),
      EchelleZ( par.echelleZ() )
{
    ReduitCoordonnes = MatriceGL3( 1. / EchelleX, 1. / EchelleY, 1. / EchelleZ );

    MoteurResolution = new Solveur();
}

/** Destructeur, virtuel pour assurer une séquence de destruction correcte.
 *
 * Destructeur, qui doit libérer la mémoire allouée pour le solveur
 * numérique.
 */
BaseHyperboloide::~BaseHyperboloide()
{}

/** Calcule l'intersection entre le rayon et le cône infini formé par
 * l'instance.
 *
 * Pour se simplifier la vie, on se ramène dans le système de coordonnées
 * locales, puis on fais une homothétie qui réduit tous les facteurs d'échelles
 * à 1. Le solveur résoud ensuite l'équation du 2ème degré associée.
 *
 * @param[in] ray instance de rayon avec lequel il faut trouver l'intersection.
 *
 * @retval NullePart si l'intersection n'existe pas, @retval local la
 * coordonnée du point d'intersection.
 */
Coord BaseHyperboloide::trouveIntersection(const Rayon& ray) const
{
    Coord origine( VersLocal( ray.origine() ) );
    Vecteur dir( VersLocal( ray.direction() ) );

    origine *= ReduitCoordonnes;
    dir *= ReduitCoordonnes;

    double a( dir.x() * dir.x() + dir.y() * dir.y() - dir.z() * dir.z() );
    double b( dir.x() * origine.x() + dir.y() * origine.y() -
              dir.z() * origine.z() );
    b *= 2.;
    double c( origine.x() * origine.x() + origine.y() * origine.y()
              - origine.z() * origine.z() - 1. );

    MoteurResolution -> setParams( c, b, a );

    vector< double > solutions;

    MoteurResolution -> resoud( solutions );

    vector< double >::const_iterator It, Begin( solutions.begin() ),
            End( solutions.end() );

    Coord local;

    for ( It = Begin; It != End; It++ )
    {
        if ( *It > 1e-7 )
        {
            local = origine;
            local += (*It) * dir;
            local *= ReduitCoordonnes.inverse();
            return VersLocal( local );
        }
    }

    return NullePart;
}

/** Calcule la valeur du champ de normale au point donné.
 *
 * La normale est donnée par le gradient de l'équation de la surface. Pour
 * simplifier on se ramène dans les coordonnées locales, puis on transforme le
 * résultat obtenu.
 *
 * @warning Ne vérifie par si le point est à la surface du volume.
 *
 * @param[in] coord point de la surface.
 *
 * @return la normale au point @c coord.
 */
VecteurU BaseHyperboloide::normale( const Coord& coord ) const
{
    Coord pt( VersLocal( coord ) );

    Vecteur local( 2. * pt.x() / ( EchelleX * EchelleX ),
                   2. * pt.y() / ( EchelleY * EchelleY ),
                   -2. * pt.z() / ( EchelleZ * EchelleZ ) );

    return VecteurU( VersGlobal( local ) );
}

/** Permet de vérifier si un point est à @e l'intérieur de l'instance.
 *
 * Détermine si le point est dans le volume délimité par l'équation de
 * la surface.
 *
 * @param[in] coord point à tester.
 *
 * @retval true si le point est bien dans le volume, @retval false sinon.
 */
bool BaseHyperboloide::estDans(const Coord& coord) const
{
    Coord pt( VersLocal( coord ) );

    if ( pow( pt.x() / EchelleX, 2 ) + pow( pt.y() / EchelleY, 2 ) -
         pow( pt.z() / EchelleZ, 2 ) - 1. <= -1.e-7 )
        return true;
    else
        return false;
}

/** Permet de vérifier si un point est à la surface de l'instance.
 *
 * Détermine si le point est sur la surface.
 *
 * @param[in] coord point à tester.
 *
 * @retval true si le point est bien sur la surface, @retval false sinon.
 */
bool BaseHyperboloide::estSur(const Coord& coord) const
{
    Coord pt( VersLocal( coord ) );

    if ( abs( pow( pt.x() / EchelleX, 2 ) + pow( pt.y() / EchelleY, 2 ) -
              pow( pt.z() / EchelleZ, 2 ) - 1. ) < 1.e-7 )
        return true;
    else
        return false;
}

/** Méthode affichant les paramètres de l'instance.
 *
 * Cette méthode affiche le nom de la forme et le sommet, ainsi que les axes
 * avec le facteur d'échelle qui leur est associé.
 *
 * @param[in] os référence sur le flot de sortie à utiliser.
 *
 * @return la référence sur le flot de sortie utilisé.
 */
ostream& BaseHyperboloide::affiche(ostream& os) const
{
    return os << "Forme : " << nom() << endl << "Paramètres :" << endl
              << "Sommet    : " << Sommet << endl
              << "Axe X     : " << AxeX
              << ", facteur d'échelle : " << EchelleX << endl
              << "Axe Y     : " << AxeY
              << ", facteur d'échelle : " << EchelleY << endl
              << "Axe Z     : " << AxeZ
              << ", facteur d'échelle : " << EchelleZ << endl;
}

/** Méthode calculant les coordonnées intrinsèques d'un point de la surface.
 *
 * Si le point n'est pas sur la surface, on renvoie NullePart.
 *
 * Permet d'obtenir un point en 2D d'après la position en 3D d'un
 * point de la surface. Pour avoir une bijection, on utilise un angle
 * (définit dans le plan normal à l'axe de symétrie) et une longueur
 * (la "hauteur" du point).
 *
 * @param[in] pt point de la surface.
 *
 * @retval NullePart si le point n'est pas à la surface, @retval local dont les
 * coordonnées sont @code ( theta, z, 0 ) @endcode
 */
Coord BaseHyperboloide::intrinseque(const Coord& pt) const
{
    if ( ! estSur( pt ) )
        return NullePart;

    Coord local( VersLocal( pt ) );

    Coord retour( 0., local.z(), 0. );

    local.setZ( 0. );

    VecteurU dir( local ), x( 1., 0., 0. ), y( 0., 1., 0. );

    retour.setX( acos( dir * x ) );

    if ( dir * y < -1.e-7 )
        retour.setX( 2 * ValeurPi - retour.x() );

    return retour;
}

/** Accesseur à Sommet.
 *
 * Permet d'accéder au "centre" de l'Hyperboloide.
 *
 * @return une référence constante sur Sommet.
 */
const Coord& BaseHyperboloide::sommet() const
{
    return Sommet;
}

/** Accesseur pour AxeX.
 *
 * Permet de récupérer l'axe X de l'Hyperboloide.
 *
 * @return une référence constante sur AxeX.
 */
const VecteurU& BaseHyperboloide::axeX() const
{
    return AxeX;
}

/** Accesseur pour AxeY.
 *
 * Permet de récupérer l'axe Y de l'Hyperboloide.
 *
 * @return une référence constante sur AxeY.
 */
const VecteurU& BaseHyperboloide::axeY() const
{
    return AxeY;
}

/** Accesseur pour AxeZ.
 *
 * Permet de récupérer l'axe Z de l'Hyperboloide.
 *
 * @return une référence constante sur AxeZ.
 */
const VecteurU& BaseHyperboloide::axeZ() const
{
    return AxeZ;
}

/** Accesseur pour EchelleX.
 *
 * Permet d'accéder au facteur d'échelle selon Ox.
 *
 * @return une référence constante sur EchelleX.
 */
const double& BaseHyperboloide::echelleX() const
{
    return EchelleX;
}

/** Accesseur pour EchelleY.
 *
 * Permet d'accéder au facteur d'échelle selon Oy.
 *
 * @return une référence constante sur EchelleY.
 */
const double& BaseHyperboloide::echelleY() const
{
    return EchelleY;
}

/** Accesseur pour EchelleZ.
 *
 * Permet d'accéder au facteur d'échelle selon Oz.
 *
 * @return une référence constante sur EchelleZ.
 */
const double& BaseHyperboloide::echelleZ() const
{
    return EchelleZ;
}
